;202230163021于凯泉
;导入数据，提取数据，亮度温度计算

PRO MOD021_BT,inputfile,outputfie
  COMPILE_OPT IDL2
  ; Start the application
  e = ENVI(/HEADLESS)

  ;恢复其核心的save文件
  ENVI, /RESTORE_BASE_SAVE_FILES
  ENVI_BATCH_INIT

;  inputfile = 'D:\IDL\data\output\MOD021KM.A2019106.0120.061.2019106132509_reflectance_georef.dat'
;  outputfie = 'D:\IDL\data\output\MOD021KM.A2019106.0120.061.2019106132509_BT.dat'
  
  
  envi_open_file,inputfile,r_fid=fid_nor
  if(fid_nor eq -1)then return
  ENVI_FILE_QUERY,fid_nor,ns=ns_nor,nl=nl_nor,nb=nb_nor,dims=dims_nor
  map_info=envi_get_map_info(fid=fid_nor)
  pos_nor=lindgen(nb_nor)

  ;read 1 2 7 bands
  ref1=envi_get_data(fid=fid_nor,dims=dims_nor,pos=0)
  ref2=envi_get_data(fid=fid_nor,dims=dims_nor,pos=1)
  ref7=envi_get_data(fid=fid_nor,dims=dims_nor,pos=6)

  ;read 21 31 32 bands
  rad21=envi_get_data(fid=fid_nor,dims=dims_nor,pos=23)
  rad31=envi_get_data(fid=fid_nor,dims=dims_nor,pos=32)
  rad32=envi_get_data(fid=fid_nor,dims=dims_nor,pos=33)
  
  ;brightness temporature
  wv=[3.99157,11.0121,12.0159]
  tcs=[0.9998646,0.995608,0.9997256]
  tci=[0.09262664,0.1302699,0.07181833]

  tb21=cal_tb(temporary(rad21),wv[0],tcs[0],tci[0])
  tb31=cal_tb(temporary(rad31),wv[1],tcs[1],tci[1])
  tb32=cal_tb(temporary(rad32),wv[2],tcs[2],tci[2])

  BTarr = fltarr(ns_nor,nl_nor,6)
  BTarr[*,*,0] = ref1
  BTarr[*,*,1] = ref2
  BTarr[*,*,2] = ref7
  BTarr[*,*,3] = tb21
  BTarr[*,*,4] = tb31
  BTarr[*,*,5] = tb32

  ENVI_WRITE_ENVI_FILE,BTarr,out_name=outputfie,map_info=map_info

  ENVI_FILE_MNG, ID = fid_nor, /REMOVE
  
  ; Close the ENVI session
  e.Close

  ENVI_BATCH_EXIT

END


;亮温计算函数文件
;coculate the lightness temperature of thermal
function cal_tb,radiance,wv,tcs,tci

  C1=1.1910659e8
  C2=1.438833e4
  x1= wv*alog(1+c1/(wv^5*radiance))
  result=(C2/x1 - tci ) / tcs
  return,result
end